MTPA を解くMATLABプログラムのサンプル
埋込型永久磁石同期機 (Interior Permanent Magnet Synchronous Machine: IPMSM) のトルク最大化 (Maximum Torque Per Ampere: MTPA) を数値的に解くMATLABプログラム。fminconを用いて非線形最適化問題を解くため，Optimization Toolbox が必要となる。
$ T = p\left\{\Psi_{m}\, I_{q}+L_{d}\left(1-\xi\right) I_{d}\, I_{q}\right\} \longrightarrow \max
Subject to:
$ \left(I_{d}+\frac{\Psi_{m}}{L_{d}}\right)^{2}+\xi^{2} I_{q}^{2} \leq \frac{V_{\max}^{2}}{\omega^{2} L_{d}^{2}} (Voltage limit)
$ I_{d}^{2}+I_{q}^{2} \leq I_{\max}^{2} (Current limit)
where
$ I_{d}, I_{q}: Current on d-q frame
$ T: Torque
$ p: Number of pole pairs
$ \Psi_{m}: Flux by PM
$ L_{d}: Self inductance of d windings
$ \xi = \frac{L_{q}}{L_{d}}: Saliency
$ V_{\max}: Maximum voltage
$ I_{\max}: Maximum current
$ \omega: Angular frequency
code:MATLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Program to maximize the torque of PMSM under specific current
% for Maximum Torque Per Ampere (MTPA) rule
% developed for MDS at Sophia Univ. by Prof. M. Miyatake
%
% Required toolbox (must be installed):
% "Optimization Toolbox"
%
% 2021/06/20 Ver.1.0 program & comments written
% 2023/06/19 Ver.1.1 revised and merged to 1 file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Motor specifications
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p = 5; % Number of pole pairs
Psi = 1.0; % Flux of PM Wb
Ld = 0.05; % Inductance at d axis H
xi = 3; % Saliency (xi = Lq/Ld)
Vmax = 240; % Maximum Voltage V
Imax = 20; % Maximum Current A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Motor operation constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 200; % Rotation speed min-1
omega = 2*pi*(N/60)*p; % Angular frequency rad/s
accel = true; % if accel = true, acceleration, else deceleration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameter initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ivec = 0, 0; % Initialization of Id & Iq A
Il = -Imax, -Imax; % Lower bound of Id & Iq A
Iu = 0, Imax; % Upper bound of Id & Iq A (assuming xi>=1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Torque maximization
% by fmimcon for nonlinear optimization with constraints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ivec, fval, exitflag = ...
fmincon(@(x)torque(x,p,Psi,xi,Ld,accel)*(1-2*accel),Ivec,[],[], ...
[],[],Il,Iu,@(x)viconst(x,p,Psi,xi,N,Ld,Vmax,Imax,accel,omega), ...
optimset('Algorithm', 'sqp', 'TolX', 1e-6, 'TolFun', 1e-6, ...
'MaxFunEvals', 1000, 'MaxIter', 1000 ));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Results and graph drawing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Id = Ivec(1)
Iq = Ivec(2)
Torque = torque(Ivec,p,Psi,xi,Ld,accel)
angle = (0:2:360)*2*pi/360;
di = cos(angle)*Imax;
qi = sin(angle)*Imax;
dv = cos(angle)*Vmax/omega/Ld - Psi/Ld;
qv = sin(angle)*Vmax/omega/Ld/xi;
dt = -2*Imax:(Imax/50):0.4*Imax;
qt = Torque/p./(Psi+(1-xi)*Ld*dt);
plot(di,qi);
xlim(-2*Imax 2*Imax)
ylim(-2*Imax 2*Imax)
title('PMSM current control strategy')
xlabel('Id A')
ylabel('Iq A')
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
daspect(1 1 1)
hold on
plot(dv,qv);
plot(dt,qt);
plot(Id,Iq,'-o','MarkerSize',8);
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function T = torque(x,p,Psi,xi,Ld,accel)
%Torque is calculated from Id and Iq
% Extracting each current value
Id = x(1);
Iq = x(2);
% Torque calculation
T = p*(Psi*Iq+Ld*(1-xi)*Id*Iq);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c,ceq = viconst(x,p,Psi,xi,N,Ld,Vmax,Imax,accel,omega)
%Voltage and current constraints are evaluated.
% Extracting each current value
Id = x(1);
Iq = x(2);
% Inequality constraints
% First and second terms are voltage and current constraints
c = (Id+Psi/Ld)^2 + (xi*Iq)^2 - (Vmax/omega/Ld)^2, Id^2 + Iq^2 - Imax^2;
% If the value is not larger than 0, it satisfies the constraint.
% No equality constraints
ceq = [];
end